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ISENTROPIC FLOW SOLUTIONS FOR REACTING GAS 
MIXTURES IN THERMOCHEMICAL EQUILIBRIUM 

By Ernest V. Zoby, Jane T. Kemper, 
and Casimir J. Jachimowski 
Langley Research Center 

SUMMARY 

A machine program has been developed for the calculation of local (edge of the 
boundary layer) pressure, density, temperature, enthalpy, entropy, the mole fractions of 
the chemical species, velocity, and the derivative of the velocity with respect to the pres- 
sure over a blunt body from given free-stream conditions. The local values are computed 
for arbitrary reacting ideal gas mixtures in thermochemical equilibrium with the assump- 
tion of isentropic flow and a given pressure distribution over the body. Also, the asso- 
ciated normal shock and stagnation- point parameters are computed. Excellent correla- 
tions are obtained with existing solutions for air at free-stream velocities from 2000 to 
50 000 ft/sec (0.6096 to 15.24 km/sec) and altitudes up to 300 000 ft (91.44 km). Also, 
results for the normal- shock, stagnation- point, and local conditions based on an assumed 
Martian atmosphere are presented. The program was developed for the IBM 7094 elec- 
tronic computer at Langley Research Center. The identification of the program inputs, 
a flow diagram, a listing of the species and the program, and a sample input and output 
are given in the appendixes of the paper. 

INTRODUCTION 

For the calculation of local, blunt-body aerodynamic heat transfer and shear stress 
(for example, ref. 1), it is necessary to know the local (edge of the boundary layer) con- 
ditions over the blunt body. With known stagnation-point conditions, a given pressure dis- 
tribution over the body, and the assumption of isentropic flow, the local conditions can be 
obtained from a Mollier diagram or by a curve-fit technique. However, these processes 
of obtaining the local conditions are time-consuming and usually result in a loss of 
accuracy. 

Stagnation- point solutions with corresponding normal-shock data, Mollier diagrams, 
and curve fits to the Mollier diagram are available for an air model. (See refs. 2 to 7.) 
However, for gas models other than air (for example, Mars), these aids for computing the 
local conditions are not available. 



Because of the problems of time consumption, inaccuracies, and insufficient sources 
of data, a computer program to calculate local conditions for blunt-body reentry in arbi- 
trary gas mixtures is required. The only known programs which apply to the expansion 
of reacting gas mixtures are nozzle flow programs (for example, ref. 8). These pro- 
grams require computed stagnation values and given area ratios and, consequently, are 
not readily adapted to a blunt-body reentry problem. Therefore, a program for the cal- 
culation of local, blunt-body conditions and the accompanying normal- shock and 
stagnation-point conditions from given free-stream conditions has been developed. The 
program combines the normal- shock relations and a general thermochemical equilibrium 
program for arbitrary ideal gas mixtures (ref. 9) with the assumption of isentropic flow 
and a known pressure distribution over the body. Since the equilibrium calculations of 
reference 9 are used, the blunt-body local conditions and the corresponding heating rates 
and shear stress values can now be computed for arbitrary reacting gas mixtures. 

This paper describes the operation of the program, compares existing solutions for 
an air model with the present results, presents results based on an assumed Martian 
atmosphere, and includes program inputs (appendix A), flow diagrams (appendix B), 
listing (appendix C), and a sample input and output case (appendix D). 

SYMBOLS 


H enthalpy, ergs/g 

M molecular weight of mixture 

p pressure, dynes/cm 2 

R universal gas constant, ergs/mole-°K 

T temperature, °K 

U velocity, km/sec 

x distance from stagnation point corresponding to pressure, cm 

x^ mass fraction for species i 

yj mole number for species i, moles of species i/gram of mixture 

p density, g/cc 
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S/R entropy 

e convergence criteria in stagnation point and local flow solution 

r convergence criteria in normal-shock solution (=10"®) 

Subscripts: 

a denotes enthalpy change across shock based on equation (5) 

b denotes enthalpy change across shock based on equilibrium program 

1 conditions ahead of shock 

2 conditions behind shock 

e local flow conditions 

s stagnation- point conditions 

i denotes species 

n denotes iteration 

METHOD OF CALCULATION 

The one -dimensional steady flow of a gas through a normal-shock wave is given by 
the following equations : 


PlUl = P 2 U 2 

(1) 

Pi + Pl U l 2 = P2 + P2 U 2 2 

(2) 

h i + |u 1 2 = h 2 + |u 2 2 

(3) 


These equations describe the requirements of mass, momentum, and energy conservation. 
From these equations and the equation of state 


the following equations are obtained: 





( 5 ) 


p 2 = P X <1 + 


M ^ 2 


RT 


i L 



( 6 ) 



P 2 m 2 t i 


(7) 


The equilibrium flow calculations are performed by using an equilibrium program 
developed by Allison (ref. 9), which utilizes the partition function of statistical thermody- 
namics to determine the thermodynamic parameters and a free-energy of minimization 
technique by White (ref. 10) to determine the equilibrium composition of the gas. Also, 
the necessary thermochemical data are listed in reference 9. The thermodynamic prop- 
erties obtained with the general thermochemical equilibrium program for arbitrary ideal 
gas mixtures are compared with the results of other investigations in references 11 and 
12. The program provides a relation between the enthalpy, pressure, temperature, and 
composition 


H = H(T,p, Xi ) (8) 

which is coupled with equations (5), (6), (7) and an equation for the molecular weight of 
the gas 

M =(s yX 1 (9) 


Normal- Shock Conditions 

The conditions behind the normal shock are computed with the following inputs : the 
temperature, pressure, velocity, and gas composition ahead of the shock wave, an initial 
estimate for the temperature behind the shock, and the density ratio P^P 2 across the 
shock. 

The initial estimate of Pj/P 2 is used in equations (5) and (6) to obtain AH a and 
P 2 . This value of p 2 , an initial estimate of T 2 (which is equal to T(l)), and the gas 
composition ahead of the shock are used with equation (8) in the equilibrium program to 
obtain H(l), a first estimate for H 2 . (The free-stream conditions (Tj, pj, x^) are used 
directly in the equilibrium program to compute Hj.) The iteration technique used for 
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finding T 2 is the Newton- Raphson technique. To begin the iteration the first point T(l) 
is the initial T 2 estimate with the second point arbitrarily chosen as T(2) = T(l) + 10. 
The term AH b (2) which is equal to H(2) - Hj determined from the equilibrium pro- 
gram when T 2 = T(2) is now compared to AH a (determined from eq. (5)) through the 
convergence criterion, 


AH a - AH b(n) 
AH a 


= T = 10" 5 


( 10 ) 


If the convergence criterion is not satisfied, a correction to T 2 is obtained, as noted 
previously, by the Newton-Raphson technique, 


T . _ T - f(Tn) (11) 

T n+ 1 - T n ft(Tn) (ID 

where the functional relationship is written as 

f(T n ) = AH a - (AH b ) n (12) 

and the first derivative with respect to temperature is 


f’(T n ) 


(AH b ) n - (AH b ) n _! 
^n ~ ^n-l 


(13) 


The quantity (AH b ) n refers to the value of AH b calculated from the nth estimate value 
of T 2 - The temperature is repeatedly corrected until the convergence criterion is 
satisfied. 


A second estimate for Pi ^2 I s obtained from equation (7) by using the conditions 
behind the shock that satisfied the convergence criterion. This estimate of PijP 2 I s 
used in equations (5) and (6) to obtain new estimates for AH a and p 2 . This new esti- 
mate for P 2 and the previous estimate for T 2 are used to obtain a new AH b , and 
T 2 is corrected until the convergence criterion is satisfied. The process of obtaining 
better estimates for Pj^P 2 and T 2 is continued until a given T 2 , p 2 , and Pj^P 2 

satisfy equations (5), (6), (7), and the convergence criterion. (A flow diagram for the 
computation of the normal shock conditions is given in appendix B.) This procedure 
yields values of T 2 , P 2 > H 2 , U 2 , P 2 , (S/R) 2 , and the equilibrium gas composition. 


Stagnation Conditions 

The requirements at the stagnation point are 

H s = h 2 + | U 2 2 (14) 
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and 



( 15 ) 


Briefly, the calculation procedure is as follows. Equation (14) is used to obtain 
H s , the stagnation enthalpy. By using the pressure behind the shock wave pg and an 
initial arbitrary estimate for the stagnation temperature T 2 + 10, a value of the enthalpy 
is calculated through the equilibrium program. This enthalpy is compared with H s 
through the convergence criterion defined by e . (e is initially set equal to 0.01, and 
during the iteration procedure, it is reduced to r before the convergence is satisfied.) 


H s- H 

H s 


^ e 


(16) 


Better estimates for the temperature are obtained with the relation 


AT = — ^ (17) 
H n ~ H n _i 
T n ~ T n _i 


where H n is the value of the enthalpy calculated from the nth estimate of the tempera- 
ture. This process is repeated until the convergence criterion is satisfied for a given 
temperature. This temperature and a new estimate for the stagnation pressure (P 2 + Ap), 
where initially Ap = 0.1p2 and for subsequent corrections 


(S/R) s - (S/R) n 
(S/R) n - (S/R) n _ 1 


Pn - Pn-1 


(18) 


are then used in the equilibrium program to calculate S/R. Better estimates for p are 
obtained until the convergence criterion 


(S/R) s - S/R 
(S/R) s 


Se 


(19) 


is satisfied. This estimate for the stagnation pressure and the last best estimate for the 
stagnation temperature are used in the equilibrium program to compute H. The value of 
H is compared by equation (16) with H s . If the convergence criterion is not satisfied, 
the temperature is corrected. This separate iteration on the temperature and pressure 
is repeated until, for a given T and p, both equations (16) and (19) are satisfied. (The 
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flow diagram for the computation of the stagnation- point condition is given in appendix B.) 
This temperature T s and pressure p g are used to calculate the density and the gas 
composition at the stagnation point. 


Local Conditions Over a Body 


The inviscid flow at the outer edge of the boundary layer is considered to be isen- 
tropic. The normalized pressure distribution p e |p s , a required input, and the stagnation 
pressure are used to compute the local pressure. For each pressure the method of cal- 
culation, given in flow diagram form in appendix B, consists of an iteration on the tem- 
perature until the convergence criterion of equation (19) is satisfied. The corrections to 
the temperature are given by 


(S/R) g - (S/R) n 
(S/R) n _ 1 - (S/R) n 
" Tn-1 


( 20 ) 


The initial estimate for the temperature is arbitrarily chosen to be T s - 10. When the 
convergence criterion is satisfied, the equilibrium gas composition, the density, enthalpy, 
velocity, and the derivative of the velocity with respect to the pressure are determined. 
The last term is obtained from the inviscid momentum equation for any x and is normal- 
ized in the program with the free-stream pU product as 


a 


(pu), m = - 


(pU)] 

(pUL 


( 21 ) 


This parameter is important since it can be used to determine the local velocity gradient 


(sa 


if the pressure distribution over the body is known. The local velocity gradient is 

an important parameter in aerodynamic heat-transfer and shear-stress calculations. (See 
ref. 1.) Equation (21) does not apply at the stagnation point, but the stagnation- point 
velocity gradient can be determined by 



( 22 ) 


where R e ff is the "effective" nose radius (ref. 13). 


RESULTS AND DISCUSSION 


Normal-shock and stagnation- point solutions for an air model are shown in fig- 
ures 1 to 6. Excellent agreement is shown between the results of the present program 
and those obtained in references 4 and 14. 
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Figures 7, 8, and 9 show isentropic flow solutions for the normalized density, tem- 
perature, and enthalpy, respectively. These parameters are shown as functions of an 
assumed pressure distribution for several velocities at an altitude of 150 000 ft 
(45.72 km). The faired manually computed results were obtained with the use of a Mollier 
diagram for air. (See ref. 3.) 

Figure 10 shows a typical variation of the normalized derivative of the velocity with 
respect to the pressure as a function of the pressure distribution. These results were 
computed at an altitude of 150 000 ft (45.72 km) and a velocity of 30 000 ft/sec 
(9.144 km/sec). As noted, this term with the pressure gradient over the body can be used 
to compute the local velocity gradient. 

Figure 11 shows a typical variation of the gas species in mole fractions based on 
isentropic flow with an assumed pressure distribution. The results were computed for 
an air model (double-ionized species being neglected) at an altitude of 150 000 ft 
(45.72 km) and a velocity of 30 000 ft/sec (9.144 km/sec). Also, some typical, normal- 
shock gas compositions computed for air with the present program are compared with 
the results of references 4 and 14 in table I at altitudes of 50 000 ft (15.24 km) and 
250 000 ft (76.2 km). 

Figure 12 shows isentropic flow solutions for the normalized density, temperature, 
and enthalpy as functions of a pressure distribution based on an assumed Martain atmo- 
sphere jvM-7^X£Q2 = 0.282; = 0.718), ref. l£j. The data were computed at a veloc- 

ity of 15 200 ft/sec (4.64 km/sec) at an altitude of 300 000 ft (91.44 km) where 
Tj = 200° K and p^ = 9.39 dynes/ cm2. Also, the associated normal shock and 
stagnation- point values are given in the figure. The gas species in mole fractions for the 
conditions given in figure 12 are presented in figure 13 as a function of a normalized 
pressure distribution. 


CONCLUDING REMARKS 

A program for the calculation of local, blunt-body conditions and the accompanying 
normal-shock and stagnation- point conditions from given free-stream conditions has been 
developed. For the calculation of these conditions, the program combines the normal- 
shock relations and a general thermochemical equilibrium program with the assumption 
of isentropic flow and a known pressure distribution. Therefore, the local, blunt-body 
conditions and the corresponding heating rates and shear stress values can now be com- 
puted from known free-stream values and a pressure distribution over the body for arbi- 
trary reacting gas mixtures. 

Excellent correlations are obtained with existing solutions for an air model, and 
results for the normal-shock, stagnation- point, and local conditions based on an assumed 
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Martian atmosphere are presented. Also, the necessary program inputs and program 
listing with a sample input and output case are given in the appendixes of the paper. 

The program was developed for the IBM 7094 electronic computer. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., April 12, 1967, 

129-01-03-02-23. 
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APPENDIX A 


INPUT IDENTIFICATION 

The input is loaded by using the Fortran IV namelist. The input symbols are as 


follows : 


$NAM1 


ISPEC 

ordered list of integers selected from preceding list, for species 
desired behind normal shock 

N 

number of species 

JMOL 

ordered list of integers to correspond to components included in 
specie list 

M 

number of integers in JMOL 

YSTO 

non- zero mole numbers, y^, as guess for composition behind shock 

AMC 

cold molecular weight of gas 

PR 

pressure ratios (P local /P stagna ti 0 n) for ^ expansion 

NPR 

number of pressure ratios (=10) 

TAU 

convergence criteria — normally l.E-5 

$NAM2 


Pi 

free-stream pressure, atm 

T 1 

free- stream temperature, °K 

U 1 

free-stream velocity, ft/sec 

T 2 

guess temperature behind normal shock, °K 

R10R2 

guess density ratio behind normal shock, P-^jP^ 


One card of identification, columns 1 to 60, follows last data card. 

The thermodynamic data for the program are on tape, as listed in reference 9, and 
are read and selected according to the ISPEC and JMOL lists by subroutine Tape. 
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APPENDIX B 


FLOW DIAGRAMS 


The flow diagram for the normal- shock solution is 
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The flow diagram for the stagnation- point solution is 
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The flow diagram for the local flow solution is 



1 for pressure ratio ' 
j T ss - guess tem- I 
iiperature at ip (l 
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APPENDIX C 


PROGRAM LISTING 


Program 

The program listing is as follows : 

SIBFTC P 1 248 LIST 

DIMENSION DHB (2 ) *T(2 ) * RHO (2 ) 

D I MENS I ONSR ( 2 ) *P(2 ) *XSAVE (35* 10 ) * TSAVE ( 10 ) *PSAVE( 10 ) * 

1 HS AVE ( 10 ) « DSA VE ( 1 0 ) * ID1 (10) 

DIMENSION X(35) * XS ( 35 ) * H ( 2 ) 

DIMENSION Y ( 35 ) 

DIMENSION USAVE ( 1 0 ) *SSAVE ( 10 ) 

DIMENSION DUDPC 10 ) 

C 

INPUT 

COMMON I SPEC (35 ) * JMOL (10) * N * M * AMC * T 1 ,P1 *U1 * 

1 YSTO (35 ) *T2 *Rl 0R2,TAU*PR (10) *NPR*EPS* ICC 
COMMON/BLOCK/ I CODE ( 35 ) * F ( 35 ) * CAPM ( 35 ) *DHF0 (35 ) *L ( 35 ) * G ( 30 * 35 ) * 

1 SMLE (30 *35 ) * CAPLAM ( 30 * 35 ) * OMEG ( 5 * 30 * 35 ) ♦ A ( 1 0* 35 ) • CONR * CONPRF * 
2CONNO * CONH * CONK ♦ P I *EPS1 *NIT,EPS2* IC1 
NAMELIST/NAM1/I SPEC* JMOLiNiM, AMC* YSTO* TAU * PR* NPR , I CC/NAM2/T 1 * 

1P1 *U1 * T2 * R 1 0R2 « ICC 
ICC = 0 

DELP= 1 0000 • 

READ ( 5 * NAM 1 ) 

C 

C SUBROUTINE TAPE SELECTS FROM TAPE THE DATA FOR EACH OF THE SPECIES 

C CONSIDERED 

C 

CALL TAPE (N» I SPEC t M * JMOL ) 

WRITE (6*100) 

100 FORMAT ( 50H 1 PROGRAM FOR COMPUTATION OF LOCAL FLOW PROPERTIES 
1 1 4X ♦ 9HP • N • 1248*10X*5H(LRC)// 

230H JANE KEMPER FOR VINCENT ZOBY ) 

1 READ(5*NAM2) 

5 READ ( 5 ♦ 1 13) ID1 

1 13 FORMAT ( 10A6) 

WRITE(6*112) I D 1 

112 FORMAT! 17H1 FREESTREAM DATA///1 H « 1 0A6 ) 

WRITE (6* 1 14 )U1 *T1 *Pl 

114 FORMAT (///12H VELOCITY =E 1 5 • 8 * 5X * 1 3HTEMPERATURE =E15*8*5X* 
110HPRESSURE =E 1 5 • 8 * 5X * 9HDENS I T Y =E15.8> 

DELT = 1 0 • 

NCOUNT = 0 
C 

DO 9 1 = 1 ,N 
9 Y( I ) = YSTO ( I ) 

C COMPUTE HI 

PI =P1*1 •01325E6 
U1 =Ul*3048.E-2 
DT2=1 • 

C 

ICELL=4 

CALL E C OM (T1 *P1 *Y*H1*X* ICELLtRHOl *SR( 1 ) > 
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C APM I =0 • 

DO 10 1 = 1 ,N 

10 C APM I =CApM I + Y ( I ) 

C APM 1=1 ./CAPMI 

11 T ( 1 )=T2 

20 P2= r P 1 * ( 1 ,+CAPMI/ (C0NR*T1 )*U1**2* ( 1 *-RlOR2 ) ) 

21 ICELL=0. 

CALL ECOM CT ( 1 > * P2 * Y * H2 • X • I CELL * RHO ( 1 ) * SR ( 1 ) ) 

I F ( ICELL)29,29*28 

28 WR I TE ( 6 * 302 ) 

302 FORMAT (//38H X(I) WOULD NOT CONVERGE FOR THIS CASE) 

GO TO 165 

29 DHA= # 5*U1**2*< 1 •-R!OR2**2 ) 

DHB ( 1 ) =H2-H1 

T ( 2 ) =T ( 1 ) +DT 2 
285 ICELL=0 

CALL ECOM (T (2 ) *P2*Y*H2*X* I CELL* RHO (2 ) * SR < 2 ) ) 

CAPM2=0 • 

DO 30 I = 1 * N 

30 CAPM2=CAPM2+Y( I ) 

CAPM2=l ./CAPM2 

I F ( I CELL ) 3 1 *31*28 

31 DHB(2)=H2-H1 

TE ST =— ( DHB ( 2 ) —DHA )/DHA 

DT2= (DHA-DHB (2 ) )/( < DHB < 2 ) -DHB < 1 ) )/ ( T ( 2 ) -T ( 1 ) ) ) 

32 IF (ABS(TEST )-TAU >40*40 *35 

35 T ( 1 ) =T ( 2 ) 

DHB ( 1 ) =DHB (2 ) 

T ( 2 ) = T ( 1 ) +DT2 

IF CT (2 ) >36*36*37 

36 T2=T2- (T2-T1 )/2* 

GO TO 1 1 

37 IF (NIT-NCOUNT >300 *300*51 
51 NCOUNT = NCOUNT + 1 

GO TO 285 

40 R1 0R2=P1/P2*CAPMI /CAPM2*T (2 )/Tl 

P2=Pl* ( 1 .+CAPM 1/ (C0NR*T1 )*U1**2* C 1 .-R10R2 > ) 

I CELL=0 

CALL ECOM (T (2) *P2*Y*H2*X* ICELL*RH0(2) *S0R2) 

IF ( ICELL*NE«0) GO TO 28 

CAPM2=0 

DO 41 1 = 1 ,N 

41 CAPM2=CAPM2+Y( I > 

CAPM2=1 ./CAPM2 
DHA=.5*U1**2*( 1 • — R 1 0R2**2 > 

DHBC2 >=H2-H1 

DT2= 1 0 • 

TEST = - < DHB C 2 > — DHA )/DHA 
1 F ( ABS ( TEST AU ) 60 * 66 ♦ 35 

300 WR I TE ( 6 * 30 1 )NIT 

301 FORMAT C//32H THIS CASE NON-CONVERGENT AFTER I4*UH ITERATIONS) 
60 U2=RlOR2*Ul 

HS=.5*U2**2+H2 
P2=P2/1 • 0 1 325E6 
WR I TE ( 6 • 220 ) 
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220 FORMAT (25H0 NORMAL SHOCK PROPERTIES) 

WR I TE ( 6 * 1 1 4 ) U2 * T ( 2 ) , P2 , RHO ( 2 ) 

WRITE ( 6 * 223 ) ( I CODE ( I),X(I),I=1,N> 

223 F ORMAT (//4X» 35HCOMPOS I T I ON OF GAS (MOLE FRACTIONS)// 

1 ( 1 IX, A3,3X*E15.8 ) ) 

WR I TE < 6 , 222 > H2 , S0R2 

222 FORMAT (//2X, 10HENTHALPY =E1 5 . 8 , 1 7X , 9HENTR0P Y =E15.8) 

WRITE (6,208) 

LET HS = STAGNATION ENTHALPY AND SORS = STAGNATION ENTROPY 
I F ( ICC-1 ) 75,165,75 
75 CONTINUE 

S0RS=S0R2 

LET TEMPERATURE BEHIND SHOCK BE FIRST ESTIMATE 
OF STAGNATION 
T A 1 =T AU 
EPS=.01 
T ( 1 ) = T ( 2 ) 

H ( 1 )=H2 

P( 1 )=P2*1 .01 325E6 
PS = P ( 1 ) 

DELP=PS*. 1 
SR ( 1 ) = S0R2 

T ( 2 ) = T ( 1 ) +DELT 

110 CALL ECOM (T (2) ,P( 1 ),Y,H(2),X, I CELL , RHO * SR ( 2 ) ) 

IF ( ICELL ) 1 05, 1 05,28 

105 SLOPE = ( H ( 2 ) — H ( 1 ) )/(T(2)-T(l ) ) 

IF ( ABS ( (HS-H (2 ) )/HS )-EPS ) 115,115,111 

111 T ( 1 ) =T (2 ) 

H ( 1 ) =H ( 2 ) 

SR ( 1 ) = SR ( 2 ) 

T (2 )=T ( 1 )+ (HS-H (2 ) ) /SLOPE 
GO TO no 

TS IS FINAL TEMP. ON CONSTANT PRESSURE CURVE. 

THIS NOW BECOMES CONSTANT FOR PRESSURE 
ITERATION. 

1 15 TS = T (2 ) 

BEGIN CONSTANT TEMPERATURE ITERATION 

116 H ( 1 ) =H ( 2 ) 

SR ( 1 ) =SR ( 2 ) 

P( 1 )=PS 

P ( 2 ) =P ( 1 )+DELP 

117 CALL ECOM (TS,P(2) • Y , H ( 2 ) ♦ X , I CELL , RHO , SR ( 2 ) ) 

IF ( ICELL ) 175, 1 75,28 
1 75 SLOPE = (SR (2 ) -SR ( 1 ) ) / ( P ( 2 ) -P ( 1 ) ) 
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IF ( ABS ( (SORS-SR (2 ) )/SORS >-EPS ) 120* 120* 1 18 
118 DP= CSORS-SR (2 ) l/SLOPE 
H ( 1 ) «H ( 2 ) 

SR ( 1 ) =SR ( 2 ) 

P( 1 )=P(2 ) 

P ( 2 ) =P ( 1 ) +DP 
GO TO 117 
C 

C PS IS FINAL PRESSURE ON CONSTANT TEMPERATURE 

C CURVE. THIS NOW BECOMES CONSTANT FOR 

C TEMPERATURE ITERATION IF SOLUTION HAS NOT 

C BEEN FOUND. 

C 

120 PS = P ( 2 ) 

H ( 1 ) »H ( 2 ) 

SR ( 1 ) =SR ( 2 ) 

I F (DELT.GT. . 1 ) DELT = DEL T/ 1 0 • 

IF (EPS.GT.TA1 ) EPS = EPS/10. 

C 

IF (ABS ( ( HS—H (2 ) )/HS )-TAU ) 126*126,121 

C 

C BEGIN CONSTANT PRESSURE ITERATION 

C 

121 CONTINUE 
T ( 1 )=TS 

T ( 2 )=T ( 1 ) +DELT 
C 

122 CALL ECOM (T (2 > *PS* Y *H (2 ) *X* ICELL *RHO,SR(2>) 
C 

SLOPE= ( H (2 >-H ( 1 > >/ ( T ( 2 >-T ( 1 > > 

C 

IF (ABS ( ( HS““H (2 ) )/HS )-EPS >125,125*123 

123 DT = ( HS — H ( 2 ) ) /SLOPE 
SR ( 1 )=SR (2 > 

H ( 1 ) =H ( 2 ) 

T ( 1 > =T ( 2 ) 

T ( 2 > =T ( 1 )+DT 
GO TO 122 
C 

125 TS = T (2 ) 

IF (DELP.GT. 1 0. ) DELP= DELP/ 10. 

IF (EPS.GT.TA1 ) EPS=EPS/10. 

C 

IF (ABS ( (SORS-SR (2 ) )/SORS)-TAu ) 126*126*116 

SOLUTION FOUND 

126 PS=PS/CONPRF 
TS = TS 

130 HS = H (2 ) 

DS^RHO 
SORS=SR ( 2 ) 
us=o. 

DDS=0 . 

DO 132 1*1 iN 

132 XS ( I >=X ( I > 

C 

TSS=TS 
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IF(NPR.EQ.O) GO TO 155 
BODY EXPANSION 

DO 150 1 = 1 * NPR 

P=PR ( I )*PS*CONPRF 

LET FIRST ESTIMATE OF TEMPERATURE BE STAGNATION TEMPERATURE 
T ( 1 ) =TSS 

CALL ECOM ( T ( 1 ) »Pi Y,H( 1 ) iX» I CELL * RHO * SR ( 1 ) ) 

T ( 2 ) =T C 1 )-10. 

133 CALL ECOM (T(2) »Pi Y»H( 1 ) * X • I CELL % RHO • SR ( 2 ) ) 

SLOPE = ( SR ( 2 ) — SR ( 1 ) )/<T<2)-T(l > ) 

IF ( ABS ( (SORS-SR (2 ) ) /SORS ) -EPS ) 1 35* 1 35 ♦ 1 34 

134 DT = ( SORS-SR C 2 ) ) /SLOPE 
T ( 1 ) =T ( 2 ) 

SR ( 1 )=SR (2 ) 

T ( 2 ) =T ( 1 )+DT 
GO TO 133 

135 PSAVE( I )=P/1 .01325E6 
TSS = T (2 ) 

TS AVE ( I ) =T ( 2 ) 

HS AVE ( I )=H( 1 ) 

DS AVE ( I )=RHO 

IF ( <H< 1 >-HS) *GT.l .E-7 ) GO TO 1355 
USA VE ( I )=SQRT(2.*(HS-H< 1 ) ) ) 

UE=US A VE ( I ) 

DUDP ( I )=RHOl*Ul*<-l • / ( RHO*UE ) ) 

1356 SS AVE ( I )=SR(2) 

DO 136 J=1 *N 

136 XS A VE ( J * I ) =X ( J ) 

GO TO 150 

1355 US AVE ( I ) =0 • 

GO TO 1356 
150 CONTINUE 

WRITE STAGNATION POINT BODY DATA 

IF (NPR-6 >155*155*160 
155 WR I TE ( 6 * 200 ) 

200 FORMAT ( 1H1 *33X*40HSTAGNATI0N POINT AND BODY EXPANSION DATA// ) 

WR I TE ( 6 * 20 1 ) ( PR ( I ) * I = 1 * NPR ) 

201 FORMAT (5H P/PS * 15X*5Hi • 000 * 7X t 5 ( 4X * F7 • 3 * 6X ) »4XtF7.3) 

WRITE (6 *202) PS * (PSAVE( l )* 1=1 *NPR) 

202 FORMAT ( / 1 5H PRESSURE ( ATM ) * 6 < E 1 5 • 8 * 2X ) * E 1 5 . 8 ) 

WRITE (6 *203) TS * <T5AVE< I )* 1=1 ,NPR) 

203 FORMAT ( 1 5H TEMPERATURE * 6 < E 1 5 .8 * 2X ) * E 1 5 .8 ) 

WR ITE ( 6*204 ) HS « ( HS AVE ( I ) * I = 1 * NPR ) 

204 F ORM AT ( 1 5H ENTHALPY * 6 < E 1 5 • 8 * 2X > * E 1 5 • 8 ) 

WRITE (6 *2222 ) SORS* <SSAVE( I ) * 1=1 *NPR) 

2222 F ORMAT < 1 5H ENTROPY * 6 ( E 1 5 . 8 * 2X ) * E 1 5 . 8 ) 

WRITE (6*205) DSt (DSAVE< I )* 1=1 *NPR> 
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205 FORMAT ( 1 5H DENSITY . 6 ( E 1 5 . 8 . 2X > ♦ E 15 .8 > 

WR I TE ( 6 . 2 1 0 ) US . ( US AVE ( I ) * I = 1 , NPR ) 

210 FORMAT ( 1 5H VELOCITY 6 ( E 1 5 . 8 , 2X > . E 1 5 . 8 ) 

WR I TE ( 6 . 2 1 1 ) DDS . ( DUDP ( I > . I = 1 , NPR ) 

211 FORMAT ( 1 5H DU/DP 6 ( E 1 5 . 8 , 2X ) . E 1 5. 8 ) 

WR I TE ( 6 » 206 ) 

206 FORMAT (/33X.32HGAS COMPOSITION (MOLE FRACTIONS)/) 

DO 157 J=1 .N 

157 WRITE (6 *207) I CODEC J) , XS(J> . (XSAVECJ. I ) . 1 = 1 .NPR) 

207 FORMAT ( 6X *A3«6X*E15«8*2X»E15.8.2X*E15«8.2X»E15»8.2X.E15»8.2X* 
1E15.8.2X.E15.8) 

GO TO 165 
160 WR ITE(6*200) 

WR I TE ( 6 . 20 1 ) (PR ( I ) . I = 1 .6 ) 

WRITE (6. 202) PS. CPSAVEH ) . 1 = 1 .6) 

WR I TE ( 6 .203 ) TS . ( TSAvE ( I ) . I = 1 . 6 ) 

WR I TE ( 6 . 204 ) HS. ( HSAVE ( I ) . I = 1 .6 ) 

WRITE (6 .2222 ) SORS. (SSAVEt I ) , 1=1 ,6) 

WR I TE ( 6 ♦ 205 ) DS . ( DSAVE ( I ) » I = 1 . 6 > 

WR ITE(6.210) US. (USAVE( I ) . I =1 .6 > 

WR I TE ( 6 . 2 1 1 ) DDS . ( DUDP (I). 1=1. 6.) 

WR I TE ( 6 . 206 > 

DO 158 J=1 iN 

158 WR I TE ( 6 . 207 ) I CODE ( J ) , XS ( J ) . ( XS AVE ( J . I ) » I = 1 . 6 > 

C 

WRITE (6. 209 ) 

209 FORMAT ( 1 HI . 33X .48HSTAGNATI ON POINT AND BODY EXPANSION DATA (CONT.) 
1 ) 

WR I TE (6.215) ( PR (I ) ♦ 1=7. NPR ) 

215 FORMAT ( 5H PS/P . 1 OX . 6 ( 4X . F7 . 3 .6X ) . 4X ♦ F7 . 3 ) 

WR I TE ( 6 . 202 ) ( PS AVE ( I ) . I = 7 . NPR ) 

WR I TE ( 6 « 203 ) ( TS AVE ( I ) . I =7 . NPR ) 

WR I TE ( 6 . 204 ) ( HS AVE ( I ) . I =7 . NPR ) 

WRITE (6.2222 ) (SSAVE( I ) . I =7. NPR) 

WR I TE ( 6 . 205 ) ( DS AVE ( I ) . I =7 . NPR ) 

WRITE (6 .210) (USAVE ( I ) . I =7. NPR ) 

WR I TE (6.211 ) ( DUDP (I ) .1=7. NPR ) 

WR I TE ( 6 .206 ) 

DO 63 J= 1 .N 

63 WRITE (6. 20 7) l CODEC J) , (XSAVE (J. I ) . I=7.NPR) 

C 

165 WR I TE ( 6 . 208 ) 

208 FORMAT ( // 1 32 H 

1 

2 , 

P1=P1/1 • 0 1 325E6 
U1 =Ul/30»48 
GO TO 1 
END 


22 


APPENDIX C 


SIBFTC ECOM LIST 

SUBROUTINE ECOM (T,P,Y,H,X, I CELL ,RHO , SOR ) 

C 

C SUBROUTINE WHICH, GIVEN A TEMPERATURE AND PRESSURE, COMPUTES 

C THE THERMODYNAMIC EQUILIBRIUM PROPERTIES OF A GAS DESCRIBED BY 

C THE INPUT, 

C 

DIMENSION SMALE <30* 35 ) , X < 35 ) 

DIMENSION E (35 ) ,Q (35 ) , CAPF I (35)«R(10«10)'B(10), 

3TEMPS (10) , BSUM (11,1) , ABLOCK (11,11) ,PTEMP(35 ) , ZETA (35) , 

4ZETAPR ( 35 ) »ALAM(35) , 

5 I P I VOT ( 1 1 ) , I NDEX ( 1 1 , 2 ) , DQ I NT ( 35 ) , Q I NT ( 30 , 35 ) 

DIMENSION Y ( 35 ) 

C 

COMMON/BLOCK/ 1 CODE (35 ) ,F(35 ) ,CAPM (35 ) • DHFO (35 ) ,L (35 ) ,G ( 30, 35 ) , 

1 SMLE (30,35 ) , CAPLAM ( 30 , 35 ) , OMEG < 5 , 30 * 35 ) , A ( 1 0 , 35 ) , CONR , CONPRF , 
2CONNO , CONH , CONK , P I ,EPS1 ,NIT,EPS2, IC1 
COMMON I SPEC (35 ) , JMOL (10) ,N «M, AMC,Tl ,P1 ,U1 , 

1 YSTO (35 ) , T2 , R 1 OR2 , T AU , PR < 10 ) ,NPR«EPS, ICC 
EQUIVALENCE (SMLE (1,1) , SMALE ( 1 , 1 ) ) , ( I CODE ( 1 ) ,CODE ( 1 ) ) 

POP=P/CONPRF 

C 

PI=3#14159 
C=2*99793E10 
NCOUNT = 0 
LTEST=LTEST 
N2 = N 

34 TK=CONK*T 
RT=CONR*T 

DO 999 J=1 ,M 

B ( J ) =0 • 0 

DO 999 1 = 1 ,N 

B( J)=B( J)+A( J, I >*Y( I ) 

999 CONTINUE 

346 YB AR = 0 • 0 

DO 347 I = 1 , N 

347 YBAR = YBAR-f Y ( I ) 

DO 40 1=1 ,N 

TEMPI =0 
LEND=L ( I ) 

DO 37 L 1 = 1 , LEND 
IF (F ( I ) )31 ,35, 31 

31 PROD = 1 • 

DO 33 IC=1 , I Cl 

IF (OMEG ( I C , L 1 , I ) ) 32 , 33 , 32 

32 PROD=PROD*( 1 ,-EXP ( ~CONH*C*OMEG ( I C , L 1 , I )/TK) ) 

33 CONTINUE 
FF = F ( I ) 

PART = (T/(CAPLAM (LI , I ) *PROD ) )**FF 
GO TO 36 

35 PART = 1 • 

36 Q I NT (LI , I )=PART*G(L1 • I ) *EXP ( ~C ONH*C*SMALE ( L 1 , I >/TK) 

37 TEMPI =TEMP1+QINT (LI , I ) 

Q ( I ) = <SQRT(2.*PI/CONH*TK/(CONH*CONNO )*CAPM( I ) ) **3 ) *TK/CONPRF*TEMP 1 
IF ( Y ( I )/YBAR ) 38 ,38, 39 

38 CAPF I ( I ) =0 
GO TO 40 

39 CAPF I ( I )=Y( I ) * ( ALOG ( POP )+ALOG(Y( I ) /YBAR ) — ALOG ( Q ( I ) )+DHF0( I ) 

1/RT) 
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40 CONTINUE 

I F ( I CELL-4 ) 396 , 95 , 396 
396 DO 50 J=1 * M 
DO 50 K = 1 , M 
R ( K » J ) =0 • 0 
DO 50 1=1 *N 

50 R(K, J)=R(K, J)+A (J, I )*A (K, I >*Y< I ) 

C 

C SET UP MATRIX FOR SOLUTION OF EQUATIONS 

C 

DO 60 J=1 , M 
TEMPS ( J ) = 0 • 0 
DO 55 1=1 ,N 

55 TEMPS ( J ) = TEMPS ( J ) +A ( J • I )*CAPFI ( I ) 

BSUM ( J « 1 ) =B ( J ) + TEMPS ( J ) 

C 

C CONSTANT TERMS IN BSUM BLOCK 

C 

DO 56 K = 1 * M 
K 1 = K+1 

56 ABLOCK(J,Kl )=R(K,J) 

C 

C PI TERMS IN ABLOCK IN COLUMNS 2 THROUGH N+ 1 

C 

60 ABLOCK ( J * 1 )=B( J) 

C 

C (X/Y) TERMS IN FIRST COLUMN 

C 

Ml =M+1 

ABLOCK ( M 1 , 1 )=0.0 
DO 61 K=1 ,M1 
K 1 =K+ 1 

61 ABLOCK (Ml ,K1 )=B(K) 

BSUM (Ml t 1 )=0.0 

DO 62 I = 1 , N 

62 BSUM (Ml , 1 )=BSUM (Ml,l )+CAPFI ( I ) 

C 

C MAT I NV EXPECTS AN M+l BY M+l MATRIX 

C 

CALL MAT I NV (ABLOCK ( 1 , 1 ) ,M1 * BSUM ( 1 , 1 ) , 1 ♦ DETERM , IP I VOT ♦ INDEX ♦ 1 1 *0 ) 
C 

C RETURN WITH ANSWERS IN BSUM 

C 

ZETAP=BSUM ( 1 , 1 >*YBAR 
ZER0=0 • 

NEG = 0 *0 
DO 70 1 = 1 ,N 
PTEMP ( I ) =0 • 0 
DO 65 J=1 * M 
J1 = J+l 

65 PTEMP ( I ) =PTEMP ( I ) +BSUM ( J 1 « 1 )*A < J « I ) *Y ( I ) 

ZETA ( I )=-CAPFI ( I )+Y( I )*BSUM (1,1 )+PTEMP( I ) 

C 

C TEST FOR NEGATIVE OR ZERO ZETA 

C 

68 IF (ZETA ( I ) 169,695,70 

69 PIECE = — Y( I )/( ZETA ( I ) — Y ( I ) ) 

I F ( P I ECE ) 69 1 ,692,691 
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691 NEG=NEG+ 1 

ALAM (NEG ) =P I ECE 
GO TO 70 

692 Y< I ) = 0 
ZERO=l • 

GO TO 70 

695 IF (Y( I ) >69*70* 69 

70 CONTINUE 
C 

C FIND GREATEST NEGATIVE ZETA-Y 

C 

I F ( ZERO ) 70 0 * 70 0 * 698 

698 IF ( NCOUNT— N I T ) 699 *100*100 

699 NCOUNT = NCOUNT+ 1 
GO TO 346 

700 IF (NEG— 1 >78*71 *73 

71 ALAMPR= • 999999* AL AM ( 1 ) 

GO TO 745 

73 ARGl = AL AM ( 1 ) 

DO 74 1=2* NEG 

72 ARG2 = AL AM ( I ) 

ARG 1 = AM INI CARG1 * ARG2 ) 

74 CONTINUE 

ALAMPR= .999999* ARG 1 
745 I I C = 0 

75 ZETAP=0 

DO 76 1=1 *N 

ZETAPR (I ) = Y ( I ) + AL AMPR* ( ZETA ( I )-Y ( I ) ) 

76 ZETAP=ZETAP+ZETAPR( I ) 

DL AM = 0 

DO 77 I = 1 * N 

IF (ZETAPR ( I )/ZE TAP >77*77*765 

765 DL AM=DLAM+ (ZETA (I )-Y( I ) )*(ALOG( POP ) — ALOG ( Q ( I ) )+DHFO (I J/RT+ALO 

1G (ZETAPR ( I J/ZETAP) ) 

77 CONTINUE 

IF (DLAMJ81 *81 *80 

80 IF ( I I C — 3 >805*81 *81 
805 I IC= I IC+1 

AL AMPR = ALA MPR* • 9 
GO TO 75 

78 AL AMPR = 1 • 

GO TO 745 

C 

C CONVERGENCE TEST FOR Y(I)S 

C 

81 IF ( ALA MPR— •70)83*815*815 

815 DO 82 1=1 *N 

I F (ZETAPR ( I ) >813*816*813 
813 REL= Y ( I >— ZETAPR ( I ) 

IF i ABS ( REL ) -EPS 1 >818*818483 
818 REL = ZET APR (I ) / Y ( I ) - 1 • 

IF ( ABS ( REL ) — EPS2 >82*82*83 

816 IF(Y( I) >817*82*817 

817 GO TO 83 

82 CONTINUE 

Y ( I ) S CONVERGE 
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DO BOO 1 = 1 * N 
800 Y( I )=ZETAPR( I ) 

GO TO 95 

NON-CONVERGENCE OF Y(I)S 

83 NCOUNT=NCOUNT+l 

I F ( NCOuNT— N IT) 84 « 100. 100 

84 DO 85 1=1 ,N 

85 Y ( I )=ZETAPR ( I ) 

REPEAT WITH NEW Y(I)S AND NO. OF ITERATIONS LESS THAN NIT 

GO TO 346 
95 DO 201 1 = 1 .N 

201 X< I ) = Y ( I >*CAPM( I ) 

YB AR=0 • 0 
CAPMI =0 
DO 2026 I = 1 .N 
YBAR = YBAR+Y( I ) 

2026 CAPM I=CAPMI+X( I )/CAPM( I ) 

CAPM 1=1 .O/CAPMI 
Z=AMC/CAPMl 

ESUM=0 

DO 2029 I = 1 .N 
QSUM=0 
DQ INT ( I) =0 
LEND=L ( I ) 

DO 2028 LI =1 .LEND 
SUM = 0 

DO 2027 IC=1 . I Cl 

HOOTK = CONH*C*OMEG ( I C. LI .1 >/T« 

IF (OMEG ( IC.L1 . I > >2000.2027.2000 
2000 SUM=SUM+H00TK/(EXP(H00TK>-1 • ) 

2027 CONTINUE 

DO INT ( 1 ) =DQINT ( I >+QINT (LI . I )* (F ( I )/T* ( 1 . +SUM ) +SMALE ( L 1 . I )*CONH*C 
1 / ( TK*T ) ) 

2028 QSUM=QSUM+QI NT (LI « I > 

EC I > = 1 ./CAPM (I )* ( 1 .5*RT+RT*T/QSUM*DQINT < I l+DHFO < I ) > 

2029 ESUM=ESUM+X< I >*E C I > 

HOZRT = CAPMI *ESUM/ (CONR*T > + l . 

H=HOZRT*CONR*T*Z/AMC 

TK=T*CONK 

FSUM=0 

DO 2040 I =1 .N 

2033 IFCYC I ) >2034.2034.2035 

2034 CAPFI ( 11=0 
GO TO 2040 

2035 CAPF I ( I ) =Y ( I >* ( ALOG (POP )+ALOG ( Y ( I ) /YBAR ) -ALOG ( Q < I ) ) +DHF0 ( I ) 

l/RT) 

2040 FSUM = FSUM+CAPF I (I ) 

SOZR=HOZRT~ CAPM I *FSUM 
SOR=SOZR*Z 
RHO=P*CAPMl/RT 
OOZ=l «0/Z 
DO 300 1 = 1 » N 
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300 X < I > =X < I > *C APM I /C APM ( I ) 
I CELL- 0 
RETURN 
100 I CELL= 1 
RETURN 
END 


S I BFTC TAPE 

SUBROUTINE TAPE <N * I SPEC * J * JMOL ) 

C 

C 

C SUBROUTINE TAPE SELECTS THERMODYNAMIC DATA FROM TAPE 

C 

COMMON/BLOCK/ I CODE ( 35 ) iF( 35 ) ,CAPM (35 ) , DHF 0 (35) »L (35) « G ( 30 « 35 ) * 
1 SMLE (30t35) «CAPLAM{ 30 * 35 ) * OMEG ( 5 * 30 * 35 ) * A ( 1 0 * 35 ) % CONR f CONPRF * 
2C0NN0 * CONH * CONK * P I *EPS1 »NIT,EPS2i IC1 
D I MENS I ON BLOCK ( 1 50 ) * LB LOCK ( 35 ) * I SPEC ( 35 ) * JMOL (10)* 

1 OBL (5*30) 

READ (9 ) (LBLOCK ( I ) , I * 1 * 35 ) 

DO 1 I C= 1 * N 
I SP= I SPEC ( IC ) 

1 I CODE ( I C ) =LBLOCK (ISP) 

C 

READ (9) (BLOCK ( I ) * I = 1 * 35 ) 

DO 2 IC=1 *N 
I SP = I SPEC ( IC) 

2 F ( IC ) =BLOC< { ISP ) 

C 

tPE AD ( 9 > { BLOCK ( I) 4 1 = 1 4 35 ) 

DO 3 I C= 1 *N 
I SP= I SPEC ( IC) 

3 C APM ( I C ) = BLOCK (ISP) 

C 

READ (9) (BLOCK ( I ) * 1 = 1 *35) 

DO 4 IC=1 * N 
I SP= I SPEC ( IC) 

4 DHFO ( I C ) = BLOCK (ISP) 

C 

READ (9 ) (LBLOCK (I ) * I = 1 * 35 ) 

DO 5 I C= 1 *N 
I SP= 1 SPEC ( IC) 

5 L ( IC ) =LBLOCK (ISP) 

c 

IC=1 

DO 6 1=1*35 

READ (9) (BLOCK ( IL) * IL=1 *30) 

I F ( I SPEC ( I C ) — I ) 6 * 55 * 6 
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— 55 DO 56 LI S 1 * 30 
56 G (LI » IC)=BLOCK (LI ) 

1C*IC+1 

6 CONTINUE 

IC = 1 

DO 7 1=1*35 

READ (9) (BLOCK ( IL) * IL=1 *30 ) 

IF ( I SPEC ( IC >-I >7*65*7 

65 DO 66 Ll = h30 

66 SMLE(LI * IC)=BL0CK (LI ) 

IC=IC+1 

7 CONTINUE 
IC = 1 

DO 121=1*35 

READ (9) (BLOCK (IL)*IL=1*30) 

IF ( I SPEC ( IC ) -I >12*13*12 
13 DO 125 L I = 1 *30 
125 C A PL AM (LI * IC)=BL0CK(LI > 

IC=IC+1 
12 CONTINUE 

I I C= 1 

DO 8 1=1*35 

READ (9 > ( (OBL (IC*IL)*IC=1*5)*IL=1*30) 

IF ( I SPEC ( I IC >- I >8*75* 8 

75 DO 76 L 1=1 *30 
DO 76 IC=1 *5 

76 OMEG( IC*LI * I IC)=OBL( IC*LI > 

I IC=I IC+1 

8 CONTINUE 

IC=! 

DO 101=1*35 

READ (9) (BLOCK ( I J) * I J=1 * 10) 

I F ( I SPEC ( I C ) - I >10*85,10 
85 DO 9 I J= 1 * J 
I JM= JMOL ( I J ) 

9 A ( I J , I C ) =BLOCK ( I JM ) 

IC=IC+1 

10 CONTINUE 

CO NR =8 • 3 1 46938E7 
CONPRF = 1 • 0 1 325E6 
C0NN0=6 • 02322E23 
C0NH=6*6251 7E-27 
CONK= 1 % 38044E— 1 6 
P I =3 • 14159 
NIT = 300 
EPS 1=1 • E— 6 
EPS2=*01 
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IF ( I SPEC ( N ) — 32 >14,15,14 

14 N 1 = N-1 

IF ( I SPEC ( N 1 ) —3 1 ) 145,146,145 

145 IC1=1 

GO TO 20 

146 I C 1 =3 

GO TO 20 

15 I C 1 =4 
20 RETURN 

END 


Comments on use of ECOM Subroutine 

The program uses a routine MATINV to solve a matrix equation, AX = B where 
A is a square coefficient matrix and B is a matrix of constant vectors. Reference to 
this routine is found in subroutine ECOM following statement 62. The calling sequence of 
this routine is shown and briefly described in order to allow replacement by a similar 
routine, if necessary. 

CALL MATINV (ABLOCK( 1,1), Ml, BSUM(1,1), 1, DETERM, IPIVOT, INDEX, 11, 0) 
ABLOCK — first location of matrix A 
Ml — location of order of A, 1 < Ml ^11 
BSUM - first location of B 

I — number of column vectors 

DETERM - gives value of determinant (not used) 

IPIVOT, INDEX — temporary storage 

II — maximum order of A 

0 — factor used in computing determinant 


At return to calling program, X is stored at BSUM. 
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SAMPLE INPUTS AND OUTPUTS 


A sample input and a sample output are given in this appendix. 

Sample Input 


SDATA 

*NAM1 

PR(l)=.99,.9,.8..7».6».5,.4,.3,.2,.l , NPR = 1 0 * 

AMC= 28.962 , 

YSTO= 1 «E-l8i 1 .E— 18« 1 • E— 1 8 , 1 .E— 18,1. E— 18* 1 .E— 18. 

3.211 E— 4 « 1 .E— 18.2. 69E— 2 , 1 . E- 1 8 , 7.24E-3 . 1 .E-18, 
1 .E-18, 1 .E-18, 1 .E-18, 

JMOL = 1 .2,4.6. 

M = 4 , 

I SPEC =1 ,2, 3, 5, 6, 8, 13, 14,17, 18,19,20,21 ,22,23, 

N= 15, 

T Au = 1 • E — 5$ 

$NAM2 

P 1 =2 • 23E— 4 , T 1 =249. ,U1 =40000. ,T2=1 2300. , R 1 0R2= . 06 1 5$ 
SAMPLE CASE for AIR AT 200 , 000FT . 


Sample Output 


FR66STREAM DATA 


SAMPLE CASE FOR AIR AT 20C»0Q0FT. 


VELOCITY - 0.40C0000CE C5 
NORMAL SHOCK PROPERTIES 

VELOCITY = 0.74991 2446 05 


TEMPERATURE = 0.24900000E 03 


TEMPERATURE = 0.12343715E 05 


PRESSURE = 0.22300000E-03 


PRESSURE = 0.43624 694E 00 


DENSITY = 


DENSITY = 0 • 514881026-05 


COMPOSITION OF GAS (MULE FRACTIONS) 


E- C. 17990338E uQ 
N 0. 4906 4 7 ICE 00 
N + 0.15244G18E 00 
0 0 . 1465 743 C£ 00 
0+ 0. 26522893 fc — 01 
0- 0 .252549406— 05 
A 0.292I6158E-02 
A+ 0. 91 71 93706-03 
N2 0.42063C 18E-04 
N2 + 0 .1072454 lfc— 04 
02 0.95 20 0257 t-C7 
02+ 0.131356416-06 
02- 0.34849522E— 11 
NO 0.324051296-05 
N0+ 0.146474756-04 


ENTHALPY = 0 . 742500596 12 


ENTROPY = 0. 6745 1369E 02 
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P/PS 

PRESSURE* ATM) 

TEMPERATURE 

ENTHALPY 

ENTROPY 

DENSITY 

VELOCITY 

DU/D P 


E- 

N 

N + 

0 

L + 
C- 
A 

A + 
N2 
N2 + 
02 
02 + 
C2- 
NO 
NU+ 


l.OOu 

0.45I4C171E uu 
0.1238T397E oi) 
0.745714V3E II 
0 .6745 1 I22t u2 
0. 5302875c t-uD 
O.CGCOJCwiut-38 
C.0CC00CGUE-3 o 


0 . 1 ti 1 C3b ouE uu 
0.4db8039Dt uu 
0. 1 53395 Job uu 
<J. I46l62i:it wu 
J.266934WE-v,l 
C.25884612E-u5 
C • 29C69 l**vE-u2 
C.92o3S7/lE-u3 
u‘ • 41 72 75 02E-J4 
C. 1080785:; <_-C4 
C .9632 17;> /L-u / 
0.1 33459out-uO 
3633575ut-ii 
’) • 32 53fc8o«L-u5 
0. l46543oot:-u4 


C.990 

u • <*4688 769E 00 
0.122713388 05 
u • 74484 55 HE 12 
u. 674512 84 E C2 
u • 5 2 572 80 3E-05 
0.41657756E 05 
-C.17613590E 01 

GAS CGMPUSITICN 
U.I8C67343E 00 
u.48539726t 00 
J.153C8627E 00 
u.i46294olE OC 
0.26638 157E-01 
0.256980738-05 
U.291153G7F-C2 
U.V2348078E-03 
u • 4 1 8472 32 8-04 
u . i 0 ? 85 2 75 8-04 
o.^6C1834CE-07 
0. 13284658E-06 
u * 5 5 859 2 9CE- 1 1 
C.32511633E-05 
l .1 4655 702E-C4 


C.90C 

C.40626154E CC 
0.12258 190E 05 
0. 7366H040E 12 
0.67451 352 F 02 
0.48436431E-05 
0. 13442121E 06 

-0 . 59 30 3 694 E 00 

(MULE FRACTIONS) 
C. 177 23863E 00 
0 . 4949792 1 E 00 
0.150199518 00 
0. 14753717E 00 
0. 26 120404 E-0 1 
0 . 2 39 886 38 E-05 
0.29548934E-02 
0.89617819 E-C3 
0. 42983886 E~04 
u. 10570905E-04 
C .9316798 2 E-07 
0. 12716973E-06 
C. 32002 132 E-ll 
C. 322 63 1 95 E-C5 
0. 14666660E-04 


0.800 

C.36112136F 00 
0.121 20495E 05 
0.72673789E 12 
C.67451361E 02 
0.437b7379E-05 
0. 19481810E 06 
-0.45283693E CO 


0.17300587E 00 
0 • 501 8649 7E 00 
0 • 1466326 IF 00 
0 • 149062 72E 00 
C . 2 548 800 9 E~0 1 
G.22017069E-05 
C. 30 079672 E~ 02 
0. 862 89 164E -03 
0 * 4443 289 86-04 
0. 10308339E-04 
0.897 60 569E-07 
0. 120488836-06 
0 • 27749 7046- 1 1 
C * 3 19 5584 8E-05 
0. 146 82487E-04 


0.700 

C.31598119E 00 
0.11967188E 05 
0 * 71 56699 7E 12 
0.674515 L4E 02 
0. 390 1081 7E— 05 
0 • 245 1 3246E 06 
-0.40377164E 00 


0 . 16822970E OC 
0 . 50964249E 00 
0 • 14259949E 00 
0. 15077723E 00 
0. 24781766 E— 01 
0. 19958095E— 05 
0.30673689E-02 
0 • 8258331 6E-03 
0.46 133440E— 04 
0 ■ 100 13599E-04 
0. 860 4189 7E— 07 
0. 1 1 333713E-06 
0.235888 00 E-ll 
0. 31605238E— 05 
0. 14702887E-04 


0.600 

0.27084102E 00 
0. 1 1793714F 05 
0.7031472 2E 12 
0.67451718E 02 
0.341 53399 E-05 
0.291 77975E 06 
-0.3874651 OE 00 


0. 16274080E 00 
0. 5 1858996E 00 
0 . 1 3795477E 00 
0.15273888E CO 
0. 2 39 790 72 E- 01 
0.17797244E-05 
0. 3 1349293E-02 
0.783 95 288E-03 
0.481 83625E-04 
0.96781273E-05 
0.8 1941 868 E— 07 
0.10560910E-06 
0. 19536032E-1 1 
0. 3 120196 7E— 05 
0.14730637E-04 


0.500 

0. 22570086E 00 
0. 1 1592938E 05 
0. 68864702E 12 
0.67451230E 02 
0. 291 77761E-05 
0. 33783994E 06 
-0.39170465E 00 


0.15626631E 00 
0. 52915583E 00 
0. 1324641 7E 00 
0. 15504163E 00 
0. 23043328E-01 
0. 155 14247E-05 
0.32135655E-02 
0.7356 l 595 E— 03 
0. 50755873E-04 
0.92895 192 E-05 
0.77364895E-07 
0.97152661E-07 
0.1 56 12857 E-ll 
0. 30736785E-05 
0. 1477 1852 E— 04 


CO 
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to 


POINT AND BODY EXPANSION DATA (CONT.) 
0.200 0.100 


PS/P 

PRESSURbt ATM) 

TEMPERATURE 

ENTHALPY 

ENTROPY 

DENSITY 

VELOCITY 

DU/D P 


E- 

N 

N + 

0 

0 + 
c- 

A 

A + 

N 2 
N 2 + 
02 
C2 + 
0 2 - 
NO 
N0 + 



0 . 4u U 


c. 

18056Ct>8 6 

uG 

0. 

11353716c 

05 

; J a 

67140461b 

LZ 

c . 

674512^88 


0. 

2405587 1E- 

-o5 

0. 

3 8551 3 but 

06 

0. 

41635235b 

uO 

0. 

, 14840694b 


0. 

, 54 1 9965 96 

0 J 

c . 

,125783196 

OG 

0. 

.1578229/6 

ou 

n 

w a 

, 219222436- 

-Cl 

0, 

, L308C5 7bh- 

-o5 

c. 

> 330 7,44 55 c- 

-GZ 

0. 

.67 8 510 566- 

-ori 

0, 

,54115896b- 

-04 

0, 

.882482676 


0. 

.7212293 /6- 

-07 

c. 

.87721356b- 

-c/ 

c, 

.118435C86- 

-11 

0. 

,30175496b- 

-05 

0 

. 148315356 

- u4 


STAGNATION 
C. 300 

13542C51E OG 
1 1G54875E C5 
e4 992 539E 12 
o7451319E 02 
18747974 £-05 
43769749E 06 
47C53647E 00 

DAS COMPOSITION 
C.13E36519E 00 
u . 55 842 5 75 E 00 
0.11 722432E CC 
0.16135591E CO 
0. 20 509776 E- 01 
C. 10452274E-05 
0.34247 10CE- 02 
0.6C823576E-C3 
0.58 £671 02E-04 
0.8246123CE-05 
0.6553785 3 E-C7 
C. 765C8837E-C7 
G.62711779E-12 
0. 29482 724E-05 
0. 14928134E-C4 


0.90280 34 IE -01 
0.106 50144E 05 
0.6210080 IE 12 
0.674511 99 E 02 
0. 13183835E-05 
0 . 4994 1350E 06 

-0.5 8643466E 00 

l MOLE FRACTIONS) 
0 . 1243991 1 E 00 
0.58130766E 00 
0. 105286 89 E OC 
0.1 6624 105 E 00 
0. 18574141 E-Ql 
0. 75541206E-06 
0. 358249C1E-02 
0.51582090E-03 
0.66561 842 E-04 
0 . 747268 90 E- 05 
0.58247693E-07 
0.6 39 10139 E-07 
0.496 1483 6E- 12 
0. 285922C5E-05 
0.15109437 E-04 


0.45140170E-01 
0.99955046E 04 
0.57499205E 12 
0.67451323E 02 
0. 721 04105E— 06 
0. 5 84 3 3 36 IE 06 
-0.91643324E 00 


0 . 10 1 1 1340E 00 
0. 61951096E 00 
0.85326370E-01 
0.1 7434409E 00 
0.1538921 2E-0 1 
0.42269858E-06 
0.38310476E— 02 
0.37628723 E-03 
0. 8 36 0 2 24 6 E— 04 
0. 62715311 E-05 
0. 47599765 E— 07 
0 . 46600355E-07 
0 . 2044 830 2E— 12 
0. 2 74 1 2 40 9 E— 05 
0. 15581785 E— 04 


0 . 
u. 
0 • 
0 « 
c. 
0 • 
— 0* 
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TABLE I.- TYPICAL COMPARISON OF GAS COMPOSITION IN MOLE FRACTIONS BEHIND NORMAL SHOCK 


CO 

tn 





Mole fraction behind normal shock for Ui 

of - 



Species 

2000 ft/sec (0.6096 km/sec) 

15 000 ft/sec (4.57 km/sec) 

30 000 ft/sec (9.145 km/sec) 

50 000 ft/sec (15.24 km/sec) 


Present method 

Reference 4 

Present method 

Reference 4 

Present method 

Reference 14 

Present method 

Reference 14 

Altitude, 50 000 ft (15.24 km) 

n 2 

7.809 x 10"* 

7.809 x 10" 1 

6.27 x 10" 1 

6.283 X 10" 1 

8.764 X 10“ 2 

9.016x10-2 

1.218 x 10- 4 

1.494 x 10' 4 

°2 

2.097 x 10' 1 

2.098 x 10" 1 

1.11 x 10' 2 

1.225 X lO’ 2 

1.316 x 10“ 4 

1.485 X 10- 4 

4.389 x lO' 6 

7.441 x lO'® 

A 

9.3 x 10' 3 

9.324 X 10' 3 

7.909 x 10“ 3 

7.949 X 10-3 

5.063 x 10' 3 

5.077 x 10-3 

1.779 X 10-3 

1.676 x lO' 3 

NO 

5.315 x 10" 13 


5.467 x lO- 2 

5.657 X 10' 2 

5.089 x lO' 3 

5.464 x lO' 3 

4.516 XIO’ 5 

5.555 x lO" 5 

N 



1.951 X lO" 2 

1.828 x 10-2 

6.689 X 10" 1 

6.654 X 10-1 

3.957 X lO'l 

3.84 X 10-1 

O 

2.061 x lO’ 32 

4.628 x 10' 32 

2.797 x 10' 1 

2.766 X 10" 1 

2.226 X 10' 1 

2.227 X 10-1 

1.182 X 10-1 

1.158 X 10-1 

e 



5.414 x 10“ 5 

5.318 x lO" 3 

5.259 x 10-3 

5.429 X 10-3 

2.418 X lO'l 

2.487 x 10' 1 

N 2 + 



2.59 x 10-® 

2.391 X 10“ 8 

3.155 x 10” 4 

3.757 X lO- 4 

1.036 x 10- 4 

1.624 x 10" 4 

0 2 + 



1.982 X 10‘ 7 

2.012 X 10-7 

3.916 X 10-6 

4.468 X 10" 6 

5.824 x 10-6 

9.276 x 10“6 

NO+ 



5.83 x 10“ 5 

5.673 x lO' 5 

6,541 x 10" 4 

7.595 X 10-4 

0.789 x 10“ 4 

1.084 x 10- 4 

N+ 



5.554 x 10-9 

4.712 X lO” 9 

3.596 x lO' 3 

3.598 X 10-3 

1.989 X 10-1 

2.05 x 10" 1 

o + 



10.308 x 10-8 

9.604 X 10' 8 

7,24 XIO" 4 

7.407 x 10- 4 

4.121 X lO’ 2 

4.203 x lO' 2 

O' 



4.315 x 10-6 

3.877 X 10' 6 

5.765 X 10' 5 

7.146 x 10-5 

2.212 x lO" 4 

4.284 x 10' 4 

A + 



5.006 x 10-10 


2.31 x 10' 5 

2.259 x 10-5 

1.765 X 10-3 

1.836 X 10" 3 

0 2 “ 



1.783 X 10-7 


6,612 x lO" 8 


2.64 x 10-8 

L 


Mole fraction behind normal shock for Uj of - 


Species 

5000 ft/sec (1.52 km/sec) 

15 000 ft /sec (4.57 km/sec) 

30 000 ft /sec (9.145 km/sec) 

50 000 ft/sec (15.24 km/sec) 


Present method 

Reference 4 

Present method 

Reference 4 

Present method 

Reference 14 

Present method 

Reference 14 

Altitude, 250 000 ft (72.20 km) 

N2 

7.808 x 10" 1 

7.808 x 10" 1 

5.879 x 10- 1 

5.898 x 10' 1 

3.184 x lO- 3 

3.386 X 10-3 

8.032 X lO’ 7 

8.865 X 10- 7 

02 

2.096 x 10- 1 

2.097 x 10- 1 

7.189 x 10-5 

8.065 X lO” 5 

2.908 x lO' 7 

3.173 x 10-7 

3.434 X lO* 9 

4.08 x lO" 9 

A 

9.301 x 10-3 

9.324 x lO' 3 

CO 

< 

0 

I— < 

X 

7.459 x 10- 3 

4.624 x 10-3 

4.626 x 10-3 

1.098 X lO’ 3 

1.126 x 10' 3 

NO 

2.577 x 10' 4 

2.5 x 10- 4 

2.549 x lO' 3 

2.658 x lO' 3 

3.187 x lO' 5 

3.424 x lO' 5 

8.754 x lO* 8 

9.884 X lO" 8 

N 

3.669 x 10~16 

3.172 x lO" 18 

7.094 x lO- 2 

6.716 x lO" 2 

7.659 x 10-1 

7.658 x 10-1 

1.94 xlO-1 

1.951 x lO" 1 

0 

7.066 x 10-7 

6.616 x 10' 7 

3.293 X 10-1 

3.328 X 10-1 

2.077 x 10-1 

2.078 x 10"! 

7.02 xlO- 2 

7.004 x lO” 2 

e 



4.378 X 10-5 

4.257 x lO- 5 

9.316 x lO' 3 

9.158 x lO" 3 

3.673 X lO' 1 

3.669 x lO" 1 

n 2 + 



9.327 x lO' 8 

8.979 x lO- 9 

0.984 x 10- 5 

1.112 x 10" 5 

1.114x10-6 

1.428 x lO' 6 

o 2 + 



4.373 x 10-9 

4.633 X lO' 9 

4.544 x 10” 8 

4.971 x 10" 8 

2.192 x 10- 8 

2.605 x lO" 8 

NO+ 



4.36 X 10-5 

4.239 x lO" 5 

6.965 x 10-5 

8.073 x lO" 5 

1.645 X 10" 6 

1.999 X lO" 6 

N+ 



1.458 X 10-8 

1.39 XIO- 8 

7.557 x 10-3 

7.379 x lO” 3 

3.024 XIO' 1 

3.016 x lO' 1 

0+ 



1.586 XIO' 7 

1.58 x 10" 7 

1.66 x lO- 3 

1.668 x lO" 3 

5.311 X lO -2 

6.341 x lO" 2 

0" 



3.865 x lO- 9 

3.329 X 10- 9 

1.008 x lO- 7 

1.001 x lO" 7 

4.087 X 10' 7 

5.429 x 10-7 

A + 



1.824 x 10-1° 


2.035 x 10' 5 

1.889 X lO" 5 

1.858 X lO' 3 

1.832 X lO- 3 

0 2 " 



1.746 X lO' 3 


2.391 x lO- 6 

1 

1.725 x lO' 8 

1 
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Figure 2 - Pressure ratio as a function of altitude. 
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Figure 4.- Enthalpy ratio as a function of altitude. 
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Figure 5.- Stagnation-point pressure ratio as a function of altitude. 
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Figure 9.- Normalized enthalpy as a function of pressure for isentropic flow. Altitude, 150 000 feet (45.72 km). 
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tr The aeronautical and space activities of the United States shall be 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
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